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Introduction 

A particulate two-phase flow CFD model was developed based 
on the FDNS code (Refs. 1,2,3) which is a pressure based 
predictor plus multi-corrector Navier-Stokes flow solver. 
Turbulence models with compressibility correction (Ref. 4) and the 
wall function models (Ref. 5) were employed as submodels. A 
finite-rate chemistry model (Refs. 6,7) was used for reacting 
flow simulation. For particulate two-phase flow simulations, a 
Eulerian-Lagrangian solution method using an efficient implicit 
particle trajectory integration scheme was developed in this 
study. Effects of particle-gas reaction and particle size change 
to agglomeration or fragmentation were not considered in this 
investigation . 

At the onset of the present study, a two-dimensional version 
of FDNS which had been modified to treat Lagrangian tracking of 
particles (FDNS-2DEL) had already been written and was 
operational. The FDNS-2DEL code was too slow for practical use, 
mainly because it had not been written in a form amenable to 
vectorization on the Cray, nor was the full three-dimensional 
form of FDNS utilized. The specific objective of this study was 
to reorder the calculations into long single arrays for automatic 
vectorization on the Cray and to implement the full three- 
dimensional version of FDNS to produce the FDNS -3 DEL code. Since 
the FDNS-2DEL code was slow, a very limited number of test cases 
had been run with it. This study was also intended to increase 
to number of cases simulated to verify and improve, as necessary, 
the particle tracking methodology coded in FDNS. 

Governing Equation 

The gas-phase governing equations of the FDNS module are the 
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Reynolds-averaged Navier-Stokes equations with the addition of 
particle drag forces and heat fluxes in the momentum equations 
and the energy equation, respectively. Due to the effect of 
large density differences between the particles and the 
surrounding gas, the drag force was considered to be the primary 
contribution to the inter-phase momentum exchange. The gas-phase 
governing equations are written as: 

J'Vapq/at) « 3[-/>U,q + H' ff G ij (dq/dZ i )]/3S i + S q 

where q - 1, u, v, w, h, k, « and a s for the continuity, 
momentum, energy, turbulence model and chemical species transport 
equations respectively. And, the transformation parameters and 
effective viscosity, are given as: 

J - f)/3(x,y,z) 

u, - (u/J) (3^,/aXj) 

Gj j = (3< { /3x k ) (3*j/3x k )/J 

A*eff ■ (*» + 

The source terms in the governing equations, S q , are given as: 

0 

-P, + V[A.. ff (Uj) x ] - (2/3)(M eff Vu) x + DX 
-Py + - (2/3)(f. eff VU) y + Dy 

-Px + Uj),] - (2/3)(/,. ff Vu) z + Dz 

Dp/Dt + h v + H p - u p Dx - v p Dy - w p Dz 

p(P r - «) 

/>(«/! 0[(C,+CjP r /O P r - c 2< ] 


where Dx, Dy and Dz represent the drag forces and n takes on 
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values between 1 and N. u p , v p and w p are the particle velocity 
components. H p is the rate of heat transfer per unit volume to 
the gas phase. stands for the viscous heat flux of the gas 
phase. P r stands for the turbulence kinetic energy production 
rate and is written as: 

P r - (M t /p) C (3U/3X, + du i /dx i ) z /2 - 2(3u k /ax k )V3] 

An equation of state, p » p/ (RT/MJ , is used to close the above 
system of equations. Turbulent Schmidt and Prandtl numbers, a q/ 
for the governing equations and other turbulence model constants 
are given taken from Refs. 4, 6 and 7. 

Finite Rate Chemistry Model 

For gas-phase chemical reaction modeling, a general system 
of chemical reactions is written in terms of the stoichiometric 
coefficients (i^j and and the i-th chemical species name (M^ 

of the j-th reaction as 
S M ( - 2 k u ' M,' 
i \ 


The net rate of change in the molar concentration of species 
i due to reactions j , X ijf is written as: 

x ij - (V^i) IK fj n (P<* i /Mu i )*' ii - 

and the species production rate, w s , (in terms of mass fraction) 
is calculated by summing over all reactions. 

u i ■ Myi sx ij 
j 
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where 

Myj m molecular weight of species i 

a, - mass fraction of species i 

p ■ fluid density 

K fj ■ forward rate of reaction j 

^ » backward rate of reaction j = K fj /K ej 

K ej * equilibrium constant 

« (l/RT) J< *' ii ''*' ij) exp^f,'*^’ - > 

f , ** Gibbs free energy of species i 
K f = A T 8 exp(-E/RT} 

Finally, the species continuity equations are written as: 
p D t a, - 7[ ( Meff /a a )Vaj] - w, 

where o a (assumed to be 0.9) represents the Schmidt number for 
turbulent diffusion. A penalty function is employed to ensure 
the basic element conservation constraints at the end of every 
time marching step. This is a crucial requirement for the 
numerical stability and accuracy of a CFD combustion model. This 
is accomplished by limiting the allowable changes in species 
concentrations, which are the solutions of the species continuity 
equations, for each time step such that the species mass 
fractions are well bounded within physical limits. The resulting 
limited changes are adjusted so that they are proportional to the 
species source terms. A similar chemistry approach and detailed 
turbulence submodels were reported previously (Ref. 8) . 

Particulate-Phase Equations 

A Euler ian-Lagrangian particle tracking method was employed 
in FDNS to provide effects of momentum and energy exchanges 
between the gas phase and the particle phase. The particle 
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trajectories are calculated using an efficient implicit time 
integration method for several groups of particle sizes by which 
the drag forces and heat fluxes are then coupled with the gas 
phase equations. The equations constitute the particle 
trajectory and temperature history are written as: 

DVj/Dt = (U, - V { )/t d 

Dhp/Dt = (T„ - T p )/t, - 6 „<f T p V (p p d p ) 


where U, 
v i 



a 


€ 


f 

p p 


* Gas Velocity 

« Particle Velocity 
= Particle Dynamic Relaxation Time 

- 4 dp/ ( 3 C d p c jU 5 - V,i) 

Particle Enthalpy 

= Particle Heat Capacity 
= Particle Temperature 
= Gas Recovery Temperature 
= Particle Thermal-Equilibrium Time 
= (p p d p )/ [12 Nu M /(Pr d p )] 

« Stefan-Boltzmann Constant 

- 4.76E-13 BTU/FT 2 -S-R 

= Particle Emissivity ** 0.20 — 0.31 

* Radiation Interchange Factor 

* Particle Diameter 
■ Particle Density 


c d and Nu stand for drag coefficient and Nusselt number for 
heat transfer which are functions of Reynolds number and relative 
Mach n umb er - Typical correlations are given in Refs. 9 and 10. 
Carlson and Hoglund's correlation (Ref. 9) is written as: 
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C d - (24/Re) (1 + 0.15 Re 0687 ) (1 + O/ 

[1 + H (3.82 + 1.28 e* 1,25t#/ ")/Re] 

Nu - (1 + 0.2295 Re 0 * 55 )/ 

[1 + 3.42 M (2 + 0.459 Re 0 - 55 ) /Re] 

where a - 0.427/M 4 - 63 + 3.0/Re°“. A more accurate but more 
complicated correlation for the drag coefficient is provided by 
Henderson (Ref. 10). That is, for Mach > 1, 

c d - 24 [Re + S (4.33 + exp (-0.247 Re/S) (3.65 - 1.53 T^/T) 

/(I + 0.353 T,/T) J]' 1 

+ exp(-0.5*M/Re 1/Z ) [0. 1M 2 + 0.2M 8 + (4.5 + 0.38a) 

/(I + a) ] + 0.6 S [1 - exp (-M/Re) ] 

where S ■ M(7/2) 1 * 2 is the molecular speed ratio, a = 0.03 Re + 
0.48 Re 1/2 . For Mach a 1.75, 

C d - [0.9 + 0.34/M 2 + 1 . 86 (M/Re) 1/2 (2 + 2/S 2 

+ 1.058 (T^Tj’^/S - 1/S 4 }] / [1 + 1-86 (M/Re) 1/2 ] 

And, for 1 < Mach <1.75, 

* . 

C d * C d M*1 + < 4 / 3 ) ( C dN=1.75 " C dH=l) 

which assumes a linear variation between M ■ 1 and M - 1.75. 

It has been shown that the Henderson drag law gives better 
motor performance predictions compared with test data. The 
applicability and possible improvement of the Nusselt number 
correlation is currently being actively researched (Ref. 11). 
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Details of the Particle Solution Method 

In the present two-phase flow model, an independent module 
was employed for the calculation of particle drag forces and heat 
flux contributions to the gas flow field. Subroutines for 
locating the particles and integrating their trajectories are 
called for each particle size group. The drag forces and heat 
fluxes are then saved for every grid point. These forces and 
fluxes are then used to evaluate the particle source terms in the 
gas-phase governing equations. In the present FDNS flow solver, 
two forms of the energy equation (i.e. static enthalpy form or 
total enthalpy form) can be selected. It has been found that 
although either fora of energy equation usually gives similar 
solutions, the static enthalpy equation provides better 
definition of the liquid rocket plume shear layers, as shown by 
extensive solutions made for the SSME. A determination of which 
fora the energy equation best simulates solid (two-phase) rocket 
motor plumes has not yet been made. 

Particle wall-boundary conditions are treated by using a 
specified fraction of the colliding particles which stick to the 
wall. Particles which stick result in a decreased particle 
velocity normal to the wall for that particle size fraction. 
Therefore, for the particle size fraction which locally collides 
with the wall, part of the particles stick and the other part is 
turned more parallel to the wall. Energy exchange is assumed to 
be due only to the particles which stick. This model of particle 
wall interaction can be improved, but new experimental test data 
must become available in order to do so. 


In the 2-D version of the FDNS flow solver, a fourth-order 
Runge-Kutta method was employed to integrate the particle 
trajectories. After a thorough test of the integration routine, 
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it was found that the explicit scheme can sometimes give diverged 
particle solutions when the source terms become large. 

Therefore, an implicit integration scheme was employed in the 
present model. For convenience, consider the X-component of the 
particle equation of motion. That is, 

dXp/dt * U p 
dUp/dt - A (U e - U p ) 

where A - l/t d 

U e = gas velocity 
Up = particle velocity 
X p = particle location 

In finite difference form the above equations can be written as: 

Xp (r*1> _ Xp (n) « (At/2) [U p <m1) + U p <n) ] 

Up (n*1) _ Up(n> „ AtA [U c - U p <n+1) ] 
or 

Xp(nH) _ Xp <n) + At/2 [U p ( "* 1} + U p (n> ] 

U (n*U . [U <n> + AtA UJ/(l+AtA) 

P P C 

These two equations are unconditionally stable despite the 
magnitude of the source terms. To provide better time 
resolution, a variable time step size is chosen so that a 
particle would take at least 4 time steps to go across a grid 
cell. 


The recognition that an improved integration scheme was 
needed for calculating the particle trajectories was a major 
hurdle in developing FDNS-3DEL. The explicit scheme appeared to 
give acceptable solutions, but detailed comparisons to previous 
FDNS-2DEL analyses showed that unacceptable pressure losses were 
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predicted. Several other factors were initially suspected of 
causing this solution behavior. Namely, the turbulence model, 
the form of the energy equation, and the particle drag law were 
initially suspected, and lengthy calculations were made before 
these effects were found not to be the cause of poor results. 
Since the FDNS-2DEL results were found to give good pressure 
field comparisons to conventional nozzle and plume flowfield 
codes (RAMP, SPP, and SPF-II) , the Runge-Kutta method was not 
expected to perform poorly in the FDNS~EL code. Resolving this 
problem consumed much of the resources which otherwise would have 
been used to run a wider variety of test cases. 
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TP.st Cases 

Th« major test case which was studied was the Tomahawk solid 
rocket motor nozzle analysis, consideration of a plume flowfield 
and of an oxygen-hydrogen coaxial injector was also made. These 
cases are described in the following paragraphs. 


• The Tomahawk Nozzle Flowf ield 

The Tomahawk nozzle flowfield was calculated with FDKS-3DEL 
and is shown in Figs. 1-4. This test case was chosen because 
comparable predictions with the FDNS-2DEL and RAMP codes had 
already been performed, and these other solutions were available 
for comparison (Ref. 12). Figures 1-4 show the velocity, Mac 
number, temperature, and water concentration profiles, 
respectively, for the chamber, nozzle, and near plume, .ne 
chamber flow was approximated to be uniform so that direct 
comparisons with the previous solutions could be made. The FDNS- 
2DEL solution predicted somewhat lower exit plane centerline gas 
temperatures (2250 «X) than the RAMP solution (2400 °K) . The 
FDNS-3DEL (2470 K) and RAMP solutions show essentially the same 
exit plane centerline gas temperatures. The raggedness in the 
temperature profile near the centerline in the nozzle appears to 
be due to a weak oblique shock. An apparent non-zero temperature 
normal gradient at the centerline in the subsonic portion of the 
nozzle flowfield is indicated. This is due to a very strong 
effect of the inlet particle flowfield boundary condition. In a 
complete SRM simulation which includes the burning grain, more 
particles would flow down the centerline from the chamber (as 
compared to the uniform flow case) and this subsonic temperature 
contour would probably change shape. The sharp breaks in the 
velocity, Mach number, and temperature contours locate the 
approximate limiting streamline of the particle laded flow with 
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respect to gas only flow which fills the nozzle. Both the static 
and total enthalpy forms of the energy equation give the same 
nozzle solutions. Letting the particles which hit the wall stick 
or elastically reflect give well behaved solutions. The only 
place where there is significant particle impact on the wall is 
at the start of the converging section. The analysis allows 
particles which hit the plume centerline to spectrally reflect, 
in order to account for particles crossing the plume centerline. 
However, particle drag moves the particles very parallel to the 
gas streamlines in the transonic region of the nozzle, such that 
such reflection does not occur in the case being considered. 

• The Tomahawk Plume Flowfield 

The near plume appears to be well predicted with FDNS-3DEL. 
The predicted free shear layer is sharply defined and indicates 
water production from afterburning reactions. Both the static 
and total form of the energy equation were considered. The total 
form of the energy equation indicates a temperature spike at the 
inception of the free shear layer. Better definition of the 
induced flow on the outside of the nozzle would probably 
eliminate such a spike. The static form of the energy equation 
does not exhibit this effect. A Mach number correction to the 
k-« turbulence model was used for this simulation. 

When the Tomahawk plume is calculated for a long distance 
down stream of the exit plane with FDNS-2DEL, excessively rapid 
plume/atmosphere mixing is predicted. This was believed to be 
due to the effect of crossing the Mach disc in the plume and 
thereby creating too much turbulent kinetic energy with the k-« 
turbulence model being used. A similar problem exists when using 
the SPF/II standard JANNAF plume code (Ref. 13). The remedy in 
the SPF/II code is to switch turbulence models between the near 
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and far plumes. An insufficient number of test cases have been 
run with FDNS-3DEL to determine if the Mach number modified 
turbulence model will indeed fix this problem, although the 
solution is better behaved than when the extended k-« turbulence 
model is used. This potential plume prediction problem for f air- 
field analysis must be left for future resolution. The FDNS-3DEL 
code should not require any change other than turbulence model 
parameters to adjust the rate of plume/ atmosphere mixing. It 
should be noted that the computed results with FDNS-2DEL, at 
first glance, look like the afterburning combustion reaction 
rates are too slow. Actually, so much of the cold atmosphere had 
mixed with the plume that the existence of afterburning was not 
apparent . 

• Liquid Injector Flowfields 


The current version of FDNS— 3DEL does not treat mass 
transport from the particle phase to the gas (or continuous) 
phase. Also, the particle phase is treated with a lumped model 
such that the particle temperature is constant throughout the 
particle at any instant of time during the flow through the 
computation field. These restrictions should be removed before 
the code is useful for describing spray combustion. However, the 
spread of a droplet cloud of supercritical fuel or LOX could be 
described with the code without modification, if one is content 
with not describing local mixture ratio changes, i.e. one assumes 
that the supercritical lump remains a lump (or particle) in the 
region of the flow being analyzed. The energy transfer for 
supercritical injection could be easily treated in this manner 
because the heat of vaporization does not have to be considered. 
In fact, models which are based on arbitrarily supplying such 
heat of vaporization (Ref. 14), do not realistically describe 
supercritical spray phenomena. The only reason that such models 
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work at all is that the heat of vaporization evaluated at the 
temperature of the oxygen lump crudely approximates the high heat 
capacity of the liquid-like lump at supercritical pressures. An 
oxygen spray eminating from a single coaxial injector could be 
described with FDNS-3DEL by assuming the oxygen lump to be of a 
constant size and density. A demonstration calculation of this 
nature was considered, but the lump density would be such a very 
strong function of the mean lump temperature that the calculation 
was not performed. If accurate real-gas equation of state models 
were used, the stated oxygen spray simulation would be 
meaningful. Currently, SECA is developing the more general 
property evaluation for a hydrogen-oxygen engine heat transfer 
analysis (Ref. 15) . However, the currently feasible constant 
property analysis was not made, because a reliable two-phase, 3- 
dimensional FDNS-3DEL code was not completed early enough in this 
study . 

Closure 

The calculation of two-phase reacting flows at best is a 
slow process. Several strategies were tried to make this process 
more efficient. Initially, ideal gas flow was computed, then the 
reactions were turned on, and finally the particle trajectories 
were calculated. The entire flowfield was calculated for each of 
these flow conditions. Recently, all of these conditions have 
been treated simultaneously from the beginning of the analysis. 
This procedure works well and results in an overall reduction in 
computation time. For analyzing rocket motors and their 
attendant plumes, it is recommended that the flowfield should be 
broken into subregions for analysis, in order to use the optimum 
step size for the Mach number range within the region. Such a 
restart option has been incorporated in FDNS-3DEL. For example, 
the motor and nozzle should be analyzed first. The computed 
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nozzle exit conditions should be used to calculate the near 
plume. The far plume should then be computed. The break between 
the near and far plume should be chosen somewhere between the 
establishment of the complex near field shock structure and the 
essentially balanced jet, predominately mixing dominated far 
field. The development of a parabolized version FDNS-3DEL to 
initially predict large plume structures and other large 
flowfields is also recommended. 
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Conclusions 

1. A two-phase, finite-rate CFD code (FDNS-3DEL) was developed 
and vectorized. The Tomahawk nozzle test case indicates the 
CFD solution accurately simulates this flow. 

2. Particle mass transfer effects are not currently included in 
the current code. The inclusion of these effects would be 
relatively simple. 

3. More test cases should be run to establish the range of 
validity of the calculation procedure. The mechanics of the 
Euler-Lagrange calculation appear to be in good working 
order. Secondary effects, such as turbulent-mixing/shock- 
structure interaction require further study with more test 
cases. However, it should be noted that suitable 
experimental data to verify many of these complex flow 
interactions are not now available. The best one can 
currently do is compare CFD solutions to SPF-II type 
analyses. 

4. Analyzing large, complex flowfields with any two-phase, 

finite-rate CFD code is a time consuming process, therefore 
utilization of all methods which would expedite such 
analyses should be considered. Analyzing the flowfields 
with carefully selected subregions and developing 
parabolized versions of the CFD codes are two such 
computational aids which should be employed. 
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